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RESUMEN 

La magnitud de las fluctuaciones de temperatura (t 2 ) que se necesitan para explicar las incoherencias observadas 
entre las metalicidades inferidas por las lfneas de recombination o por las lfneas prohibidas, no pueden alcanzarse 
con modelos al equilibrio y estacionarios en el tiempo. Por otra parte, si la fuente ionizante fuera variable, las 
fluctuaciones de temperaturas t 2 serian mucho mayores. Investigamos la respuesta temporal de la estructura 
en ionization y temperatura de una nebulosa fotoionizada por una fuenta variable periodica. Estudiamos 
como el valor medio asimptotico, (i 2 ), se comporta en funcion del periodo y amplitud de las variaciones de 
la fuente. Encontramos que las fluctuaciones de temperatura se producen solamente en la parte externa de la 
nebulosa, cerca del frente de ionization, dentro de un espesor geometrico correspondiente al 8-20% del tamano 
de la capa ionizada. Concluimos que la amplitud de variation de la estrella excitadora que se requiere para 
conseguir un (t 2 ) = 0.025 (como en la nebulosa de Orion) es excesivamente grande. La variabilidad de la 
fuente ionizante no es por lo tanto un mecanismo viable para explicar los valores observados de t 2 . Llegamos a 
conclusiones similares cuando estudiamos la variabilidad temporal que resulta de sombras intcrmitentes detras 
de condensaciones de gas opacas. Encontramos que nebulosas fotoionizadas son en promedio menos masivas 
pero algo mas calientes en el caso de fuentes variables ticlicas. 

ABSTRACT 

The magnitude of the temperature fluctuations (t 2 ) required to explain the observed inconsistencies between 
metallicities inferred from recombination lines and from forbidden lines cannot be attained by steady-state 
equilibrium photoionization models. If on the other hand the nebular ionizing source was variable, the temper- 
ature fluctuations t 2 would be significantly larger. We investigate the time-dependent response of the nebular 
ionization and temperature structure when photoionized by a periodically varying source. We study how the 
asymptotic mean value, (t 2 '), behaves as a function of the period or amplitude of the source variability. We 
find that the temperature fluctuations occur only in the outer section of the nebula, close to the ionization 
front, within a zone corresponding to 8-20% of the ionized layer's thickness. We conclude that the amplitude 
of the exciting star variations required to achieve a (t 2 ) = 0.025 (as in the Orion nebula) is unacceptably 
large. Source variability is therefore not a viable mechanism to explain the observed values of t 2 . We reach a 
similar conclusion from studies of the temporal variability resulting from intermittent shadows behind opaque 
condensations. We find that photoionized nebulae are on average less massive but somewhat hotter in the case 
of cyclicly variable ionizing sources. 
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1. INTRODUCTION 

The presence of temperature fluctuations in pho- 
toionized nebulae has been a matter of debate since 
the pioneering work of Peimbert (1967). Observa- 
tional evidence has since accumulated in favor of his 
analysis of the problematics of nebular temperatures. 
For instance, the temperatures of Hn regions and 
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planetary nebulae are observed to be significantly 
lower when derived using recombination lines rather 
than from forbidden line ratios (see review by Peim- 
bert et al. 1995). 

The work of Kingdon & Ferland (1995) and Perez 
(1997) has shown that the amplitude of the fluc- 
tuations in hydrostatic photoionization calculations 
are much smaller than required to explain the ob- 
served differences between forbidden and recombina- 
tion lines temperatures. Possible causes for the fluc- 
tuations include metallicity inhomogeneities (Torres- 
Peimbert et al. 1990; Kingdon & Ferland 1998; Liu 
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et al. 2000), photoionization by small dust grains 
(Stasinska & Szczerba 2001), and shock heating due 
to stellar wind or even supernovae in the case of giant 
Hii regions (Luridiana et al. 1999). The inclusion 
of temperature fluctuations in a theoretical frame- 
work has so far been tentative (Luridiana, Cerviho 
& Binctte 2001; Binette & Luridiana 2000; Binette, 
Luridiana & Henney 2001) given our lack of knowl- 
edge about their possible cause. 

In this Paper, we investigate how temporal vari- 
ability in ionization structure due to an intrinsically 
variable ionizing source or to intermittent shadows 
behind opaque condensations can affect the temper- 
ature structure and induce substantial temperature 
fluctuations. We also address the question whether 
a variable ionizing continuum with a fixed duty cy- 
cle would lead to a nebula that is on average warmer 
and less massive. 

2. CALCULATIONS 
2.1. Intrinsic variability of the ionizing source 

Resolution of the ionizing radiation transfer in 
the case of a variable source implies a sufficiently 
high temporal resolution in order to follow the pro- 
gression of ionizing photon fronts across the nebula. 
The atomic timescales on the other hand are much 
longer, and a lot of computation time is spent in 
the actual transfer while not much is occurring in 
terms of changes in ionization or temperature. For 
this reason, in the building of the new code yguana 
we simplified the atomic physics in order to derive 
the desired results within reasonable computation ef- 
forts. Our goal was to give precedence to a rigorous 
treatment of the transfer, and not to detailed atomic 
physics. The Courant condition, for instance, was 
rigorously satisfied in all the calculations (Eq. [201 in 
Appendix IB. 2. 1(1 . 

One of the simplifications introduced in YGUANA 
is that only Hydrogen is considered in the integra- 
tion of the nebular opacity, which is a valid approx- 
imation for Hii regions. Furthermore, we consider 
a monochromatic transfer with all ionizing photons 
having the same energy. This simplification implies 
that the hardening of the ionizing radiation with 
depth is not considered. The expected shallow ra- 
dial temperature gradient will therefore be absent 
from our calculations. This is a minor shortcoming 
of YGUANA, since the level of temperature fluctu- 
ations caused by non-equilibrium ionization largely 
exceed that produced by the temperature gradient 
alone (c.f. t |2.4fl . In effect, non-equilibrium photoion- 
ization within the IF, which we compute accurately, 
far exceeds the larger incdreased due to continuum 



hardening at large depths. The various effects in- 
troduced by UV source variability on a photoinized 
nebula will all be accurately tracked and made obvi- 
ous by comparing our time-dependent nebula to ei- 
ther the time-averaged nebula or to the steady-state 
constant-UV reference nebula. 

The cooling by metals was approximated in 
YGUANA by a fit we made to the cooling function 
using mappings ic (Ferruit et al. 1997) [c.f. Ap- 
pendix . Collisional excitation and ionization of 
Hi (and their effect on cooling) are calculated ex- 
plicitly, since they can become important within the 
ionization front (IF). The time-dependent ionization 
balance of H, the photoheating coefficients, and the 
equation of state are described in Appendix lAl while 
the practical implementation of the time-dependent 
transfer equations are presented in Appendix[5] Our 
main goal is to study the effect on the temperature 
fluctuations measured a la Peimbert as a result of 
having a variable ionizing source. Our model as- 
sumes an hydrostatic distribution of gas. 

2.2. Mean temperature Tq, t 2 and U 

Following Peimbert (1967), we define the mean 
nebular temperature, Tq, as follows 
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J v n e n H+ TdV 
J v n e n H +dV 
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where n e is the electron density, + the H n den- 
sity, T the electron temperature and V the volume 
over which the integration is carried out. The rms 
amplitude t of the temperature fluctuations is given 
by 
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J v n e n H+ (T-f ) 2 dV 
T o J v n e n H +dV 
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We have replaced tih+ by n e in the above expres- 
sions, since we do not consider explicitly He in the 
calculations, n will denote the total H density. 

We cannot compute t 2 for other ions, since we 
consider in detail only the ionization of H. We will 
drop the label H + in the following discussion unless 
required by the context. 

In the case of a variable source, both global quan- 
tities t 2 and To vary with time. The above expres- 
sions define therefore instantaneous values for a par- 
ticular time in the source temporal evolution. Even 
if the transfer of information concerning n H + , n H a 
or T from any position inside the nebula has a finite 
speed, we do not find necessary to consider this effect 
explicitly, since we are not considering any particu- 
lar external observer. Furthermore, our derivation 
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of various characteristic quantities (labeled "asymp- 
totic") will be done by integrating over a few full 
periods of the source variability. This has the ef- 
fect of washing out any phase differences introduced 
by placing a real observer at any particular external 
location. 

Let us define further useful quantities to be used 
later. One is the varying ionizing source photon lu- 
minosity, Qh- Time- averaged quantities will carry 
brackets () such as in (Qh), the time-averaged pho- 
ton luminosity. Since we adopt a spherical geometry 
and consider the source point-like, the global ioniza- 
tion parameter is defined as 

J e 2 n 2 dV f Hs £ 2„2 47rI ,2 dr V > 

where tp is the local ionizing flux (disregarding opac- 
ity), th the opacity due to H photoionization, R s the 
outer nebular radius, r the distance from the central 
source, r$ the radius of the inner cavity devoid of gas, 
e the volume filling factor, and c the speed of light. 
(This definition contemplates the possibility that e 
and the total density n may vary with radius.) 

2.3. Results for periodical sources (Model A) 

Our main goal is to derive the characteristic (t 2 ) 
of a nebula submitted to a variable ionizing source. 
The general problem is excessively vast since there 
are numerous scenarios about possible variability be- 
haviors. The perspective taken in this Paper is that 
of a Fourier analysis in which we restrict ourselves 
to the exploration of periodical sources and analyze 
the nebular response as a function of the frequency 
and amplitude of the source variability. In concrete 
terms, we have narrowed the problem to deriving the 
asymptotic (t 2 ^ for a source varying at a predeter- 
mined frequency. We define the asymptotic (< 2 ) as 
the averaged value over at least three full periods 
of the source. This has the advantage of avoiding 
a definition that depended on the initial conditions 
or on a particular ill-defined moment in the source 
history. Typically, calculations will extend to only 
five periods since we found that asymptotic values 
in most cases do not change by adding more cycles. 

A varying source generate two types of response 
in the nebular structure: a progressing IF when the 
source increases in intensity and a recombination 
front (RF) when the source decreases. These two 
phases are not symmetric. For instance, an IF has 
a finite velocity given by Vjf — tp/ns < c j where 
tp is the ionizing photon flux. Let us define y as the 
neutral fraction of H. The ionization fraction (1 — y) 



within the IF can increase at a rate, which is much 
shorter than the recombination timescale. An RF 
on the other hand propagates at the speed c, since 
a RF is initiated following propagation of the signal 
that the source is decreasing (although the neutral 
H fraction y within the front will increase only on 
timescales given by recombination). 

To illustrate these two phases, using yguana we 
will study in detail the case of source varying like a 
square-wave with a period II equal to the recombi- 
nation timescale r rec = (aBnJ^ 1 , where as is the 
recombination coefficient rate to excited states of H. 
The amplitude of the variability is ±20% (A = 0.20) 
and the duty cycle is A = 0.5. The advantage of the 
square-wave is that it will more clearly reveal the 
changes caused by either type of fronts on certain 
nebular quantities. We will assume an ionizing pho- 
ton luminosity of (Qh) = 10 51 s _1 and a nebula of 
constant density n = 200 cm~ 3 with a volume filling 
factor of e = 0.01 and an inner cavity ro surrounding 
the source of 30% of the estimated Stromgren radius, 
that is, of size ro = 4.0 x 10 19 cm. The initial physical 
conditions of the gas were given by the usual condi- 
tion of steady-state thermal and ionization equilib- 
rium, throughout the whole nebula, for a source lu- 
minosity (Qh}- The equilibrium temperature within 
the nebula turns out to be T eq = 5260 K and the re- 
combination and cooling timescales within the fully 
ionized part of the nebula are r rec — 1.1 x 10 10 s and 
Tcooi — 1.6 x 10 9 s, respectively (the cooling timescale 
is 7 times shorter than T rec ). The ionization param- 
eter (Equation [^J of the initial equilibrium model 
is U = 0.0022. This set of parameters defines our 
Model A. 

The results of the time-dependent calculations as 
a function of time r for the first five periods are 
shown in Fig. ^ The dotted line in Panel a illus- 
trates the behavior of the ionizing flux relative to 
the mean value (to be read on the upper right y- 
axis of Panel a). In all the panels, the continuous 
line corresponds to the behavior of the instantaneous 
values integrated over the whole nebula while the 
long-dashed line represents the cumulative running 
mean, which was calculated at the onset of the third 
variability cycle. We define asymptotic values as the 
last value of the running mean (after the 5th cy- 
cle is completed). The impact on the nebula of the 
progression of the IF or the recession of the RF are 
clearly visible in Panel d, which is a plot of the mass 
of ionized gas normalized to the initial equilibrium 

4 We will use primes to denote models that assume a sine- 
wave (rather than square-wave) variability, e.g. Model A" in 
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Fig. 1. All solid lines: behavior of neb- 
ular volume weighted quantities as a 
function of time r for Model A. Panel a: 
t 2 , Panel b: net cooling relative to the 
mean heating energy radiated by the 
source A net / (T heat ) [with (T heat ) = 
h(v—vo) {Qh}], Panel c: the mean elec- 
tron temperature To/T eq , and Panel d: 
the mass of ionized gas M/Al eq . In 
each panel, the long-dashed line across 
the last three cycles represents the run- 
ning mean of the plotted quantity. The 
asymptotic values referred to in the 
text (e.g. \To), (i 2 )) correspond to the 
rightmost value along the long-dashed 
line. The square-wave dotted-line in 
Panel a is the normalized photon lumi- 
nosity of the ionizing source as a func- 
tion of time, Qh(t)/ (Qh) (to be read 
on the upper right y-axis). The vari- 
ability is characterized by an amplitude 
A = 0.2, a period II = r rec and a duty 
cycle A = 0.5. 
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value. Notice that the total ionized mass of the neb- 
ula somewhat shrinks with time. Panel c shows the 
behavior of the average nebular temperature with 
time, Tq(t), normalized again to the initial equilib- 
rium value T eq . Panel b shows the behavior of the 
net integrated cooling rate of the nebula normalized 
to the the total heating available due to absorption 
of all ionizing photons (this ratio cannot exceed the 
interval ±1). With steady-state equilibrium, the net 
cooling is zero. Panel b reveals how the nebula as 
a whole heats up during the IF phase while it cools 
down during the RF phase. Finally, Panel a illus- 
trates the behavior of t 2 as a function of time. Its 
(asymptotic) average value is (t 2 ) = 0.074. 

2.4. Missing the temperature gradient 

Not considering the hardening of the UV radia- 
tion with depth results in an isothermal nebula lack- 
ing the usual outward gradient in T eq . This has 
the advantage that t 2 computed with yguana is 
determined by source variability alone. The effect 
of missing the shallow T eq gradient is negligible in 
most cases. In effect, using MAPPINGS Ic, we find 
that the constant-UV steady-state Model A is char- 
acterized 5 by a t 2 teady — 0.0026, which is generally 
very small in comparison with the values derived 
for the variable sources studied below. The true t 2 , 
which would include UV hardening, cannot be far 
off from the following estimate t 2 rue « t 2 teady + (t 2 ), 
where (i 2 ) corresponds to the time-averaged fluctua- 
tion amplitudes computed using YGUANA. We recall 
that (t 2 ) = 0.074 for Model A with A = 0.2. 

2.5. The internal nebular structure 

The changes taking place within the spherical 
nebula of Model A occur mostly within the exter- 
nal parts, that is within the outer 25%, which corre- 
sponds nonetheless to half the photoionized volume. 
This becomes apparent in Fig.|2where the H neutral 
fraction, y, and the local temperature are plotted 
as a function of nebular geometrical depth r — ro. 
The long-dashed line represents the initial equilib- 
rium model. The lack of any slope in the equilib- 
rium temperature curve (dashed-linc) in Panel b is 
the result of not considering the hardening of ioniza- 
tion radiation with depth (the transfer in YGUANA is 
monochromatic, but see H2.4fl . The 10 solid lines cor- 
respond to the ionization and thermal structure at 
equally spaced time intervals during the last (fifth) 
variability cycle (see last cycle in Fig. [TJi) . The IF 

5 The values of t 2 (0 +2 ) computed by MAPPINGS Ic in 
the constant-UV case for Models A and B are 0.0021 and 
0.0047, respectively. 



phase corresponds to the 5 curves below the dashed- 
line in Panel a while the other 5 curves correspond 
to the RF phase. During the RF phase, the neutral 
fraction lies above the dashed-line equilibrium curve 
(except near the transition to neutral gas, where 
y ~ 1). The progression of the IF is characterized, on 
the other hand, by a sharp temperature pulse above 
the equilibrium value, which is spatially resolved in 
YGUANA (see Fig.|5jD). This IF is sufficiently hot to 
cause a rise in the mean temperature (averaged over 
the whole nebula) as seen in Fig. Because the 
selected source period is rather long, being equal to 
T rec , most of the inner nebula has time to adjust its 
ionization while the source varies. This is shown by 
the coincidence of most curves for geometrical depths 
r — ?*o < 6 x 10 19 cm, as they lie (and superimpose) 
either below or above the equilibrium y dashed-line 
model in Fig. [2*. 

Interestingly, the temporal changes in the neutral 
fraction y with time for depths r — ro < 5 x 10 19 cm 
does not cause any appreciable changes in the gas 
temperature (see Panel 6). The reason is that dur- 
ing the IF or RF phase, the neutral fraction y in 
that zone evolves quickly toward the appropriate 
equilibrium (but small) value y eq over a timescale 
of order ~ 0.5y eq r rec , which is much shorter than 
T rec (because y eq < 0.01 in that zone). Hence, non- 
equilibrium heating is too short-lived in that inner 
zone for the temperature to change appreciably. 

2.6. Varying the ionization parameter U (Model B) 

The work of Campbell (1988) indicates that the 
ionization parameter 6 in H n galaxies range from 
0.0014 to 0.025. Hence Model A (U = 0.0022) is 
probably representative of nebulae 7 at the low U 
end. We have explored the behavior of (i 2 ) when 
the ionization parameter is varied. This was done 
by setting the filling factor e to unity and computing 
models which were characterized by different cavity 
sizes. By varying ro in the interval 2.5 x 10 19 < 
ro < 1.5 x 10 20 cm, we explored the following range 
0.018-0.0006 in U (we assumed a source variability 
of A = 0.05 with the other parameters such as n, A, 
n, (Qh) the same as in Model A). What we found 

6 The definition of U of Campbell (1988) is slightly differ- 
ent, which causes the values quoted by her to be 50% higher 
than those one would derive using eq. J3J 

7 Models of different densities, filling factors and source lu- 
minosities but whose product of n X e 2 X (Qh) is the same as 
well as the ratio ro / R° s , where R° s is the Stromgren radius in 
the cavity-less case (R° s = [3 (Qh) /^irnPeag] 1 ' 3 ), are homol- 
ogous and can therefore be expected to yield the same {^i 2 ^ 
values under similar variability conditions [i.e. variations of 
similar amplitude A and period II (in r rec units)]. 
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Fig. 2. Internal structure of the nebula 
as a function of geometrical depth for 
Model A (c.f. §E3J assuming a square- 
wave variability of amplitude A — 0.20. 
Panel a: the neutral fraction y of H, 
Panel 6: the local temperature T. The 
long-dashed line represents the initial 
equilibrium model. The 10 solid lines 
correspond to the ionization and ther- 
mal structure at equally spaced time in- 
tervals during the last full cycle shown 
in Fig. GJi. 
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Fig. 3. Internal structure of the nebula 
as a function of geometrical depth for 
Model B (c.f. S 12.61 assuming a square- 
wave variability of amplitude A — 0.10. 
Panel a: the neutral fraction y of H, 
Panel b: the local temperature T. The 
long-dashed line represents the initial 
equilibrium model. The 10 solid lines 
correspond to the ionization and ther- 
mal structure at equally spaced time in- 
tervals during the last full cycle. 
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Fig. 4. Behavior of the asymptotic (t 2 ^ as a function 
of the period IT for a sine-wave variability of amplitude 
A = 0.1 (Model A"). In the interval 0.2 < II < 1, (t 2 ) 
increases steeply as II 2 ' 2 . 

was that (t 2 ) increased approximately as Vu within 
the above range. 

We also studied the properties of nebulae that 
reflect more closely the geometrical particularities of 
the Orion nebulae, that is a model consisting of a 
thin sheet of ionized gas at a much higher ionization 
parameter. For instance, the Orion model of Bald- 
win et al. (1991) is characterized by a slab geometry 
at an ionization parameter of 0.03. We reproduced 
similar characteristics by using ro = 2 x 10 19 cm and 
by setting the filling factor e to unity. With (Qh) = 
10 51 s _1 and a density n — 200 cm -3 , the thickness 
of the ionized layer turns out to be 7.1 x 10 18 cm, 
that is 26% only of the final Stromgren radius. The 
ionization parameter (Equation^ of the initial equi- 
librium model is U — 0.024. This set of parameters 
defines our Model B. 

In Fig. |2| we show the internal temperature and 
ionization structure of Model B for a variability am- 
plitude of A = 0.10. This model is characterized 8 by 
an asymptotic (f 2 ) of 0.053. Overall, Models B and 
A share many similarities. For instance, the bulk of 
the variability in T or y occurs near the outer IF. 

2.7. Sensitivity to the period II 

Using a square wave variability type, we have 
illustrated the principal effects taking place within 

8 Using MAPPINGS Ic, we find that the constant-UV steady- 
state Model B is characterized by t 2 teady = 0.01 (see 112.41 . 



a nebula submitted to a periodical ionizing source. 
Two very important parameters affect the nebular 
response and the behavior of (t 2 ): the variability 
period II and the amplitude of the variations A. We 
will analyze in turns the role played by each, adopt- 
ing the same parameters as in Model A but with 
a sine- wave variability function (1 + Asin27rT/n) 
rather than a square wave. We will refer to this mod- 
ified model as Model A". We have verified that sine- 
wave models 9 with amplitude 1AA closely match the 
it 2 } from square- wave models with amplitude A. All 
periods will be expressed in units of the recombina- 
tion time, which is the natural unit for the problem 
at hand. 

Our results for Model A", assuming an ampli- 
tude A = 0.1, show a very steep dependence with 
frequency as illustrated in Fig. (t 2 ) increases 
as cx II 2 2 up to n ~ lT rec and then peaks near 
1.6r rec . The short period regime becomes progres- 
sively ill-defined in our calculations. For instance, 
the asymptotic value of (i 2 ) at II = 0.1 requires up 
to 10 cycles, because the mean values keep evolving 
well beyond 5 cycles. If radiation hardening at large 
depths had been considered, (t 2 ) would lie above the 
Steady = 0.0026 floor computed by mappings ic (see 
i J2.4H . This would alter the curve in Fig. 0] only in 
the short period domain II < 0.4. 

The main conclusion is that a nebula acts as a 
low-pass filter with only the relatively 'slow' source 
variations causing important temperature fluctua- 
tions. Clearly, n > 0.7 is the frequency domain that 
favors the largest values of (t 2 ). 

2.8. Sensitivity to the amplitude A 

Let us now vary the amplitude (up to the max- 
imum, which is A — 1 in the case of a sine- wave), 
adopting the same sine- wave variability with IL = 1. 
The results are shown in Fig. [S^,. The larger the am- 
plitude, the higher (t 2 ), as one would expect. In the 
case of Model A", there is, first a regime of faster 
increase (cx A 1 ' 8 ) of (i 2 ), followed by a slower in- 
crease regime (cx A 1 - 1 ) above A = 0.1. For Model 
B", the increase is cx A 1A everywhere. The basic 
result from this plot is that an observed t 2 bs ~ 0.02 
would require a substantial variability amplitude of 
0.09 and 0.05 from the source, for Models A" and 
B", respectively 10 , which is unacceptably high for O 
stars, as discussed in § 13.11 

9 The values of (t/ for models A" and B" in the sine- 
wave case with II = lT rec and A = 0.1 are 0.023 and 0.039, 
respectively. For A = 0.2, it is 0.055 and 0.078, respectively. 

10 However, since in Model B", t 2 steady = 0.01 (JSU, an am- 
plitude of A = 0.03 would suffice in reproducing fluctuations 
at the level (t 2 ) ~ 0.01 = 0.02 - t 2 gteady (see JTH . 
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Fig. 5. Panel a: behavior of (t 2X j as a function of the 
amplitude of the sine-wave variations of Qh- Panel b: 
behavior of (To\/T ei (curves > 1) and (M) /M eq (curves 
< 1) as a function of the amplitude of the sine-wave vari- 
ations. Model A" and B" are represented using solid and 
dashed lines, respectively. The dotted line represents the 
equation (1 + 0.7863(t 2 )) -1 ' 177 (T o )/T e9 (see § Eg) us- 
ing values from Model A". Model B" was not computed 
beyond A = 0.4. 



An interesting phenomenon occurs with variable 
ionizing sources is that the mass of ionized gas 
shrinks as (f 2 ) increases. For instance, for A = 0.40 
in Model A", the ionized mass is reduced to 91% 
of that of the initial steady-state equilibrium model. 
The reason is simple if we recall that temperature 
fluctuations (Peimbert 1967) in the case of a process 
like recombination, which is weighted towards lower 
temperatures (oc T~ - 85 ), will lead to a higher effi- 
ciency of the recombination rate, hence to a lower 
mass of ionized gas in the nebula for a given mean 
ionizing flux. (The time-averaged luminosity of all 
recombination lines, however, remains the same.) 
Note finally the tendency of the mean temperature 
(To) to increase, at large amplitudes, slightly above 
the equilibrium value, as a result of the photoheating 
energy being absorbed more efficiently during the IF 
phase. The dotted line shows the approximate be- 
havior expected for the ionized mass using the func- 
tion (M) /M eq ~ (1 + 0.5 a(a - l)(t 2 )) l ' a {f Q ) /T eq 
with a = —0.85 (adapted from Equation 5 in Binette 
& Luridiana 2000). This expression becomes invalid 
above (t 2 ) > 0.2. 

3. DISCUSSION 

3.1. Intrinsic source variability 

We have shown that variations of the ionizing lu- 
minosity can produce appreciable temperature fluc- 
tuations across the nebula. However, a careful analy- 
sis of our results leads us to believe that the proposed 
mechanism for generating fluctuations is not viable. 
For instance, to reproduce a t 2 obs — 0.025 as in the 
Orion nebula (Esteban et al. 1998) would require 
reproducing (t 2 ) = 0.015 = 0.025 - t 2 steady . This 
occurs with Model B" (more appropriate to Orion) 
when A = 0.04, that is, it requires a variability of 
±4% of the ionizing luminosity of the exciting star 
9 1 C Ori (HD 37022), assuming II ~ 1. Such vari- 
ability of Qh would imply variations of order ±250 K 
in the stellar atmosphere (derived from MAPPINGS Ic 
using T* = 39500). In the optical domain, the con- 
tinuum variations would be a lot less, about ±1% in 
amplitude (Am_g = 0.02 min-to-max). Furthermore, 
these estimates are based on a variability timescale 
very close to the recombination timescale, which is 
of order 20 years for Orion (n e ~ 5000 cm -3 ), oth- 
erwise if it was shorter it would require variation 
amplitudes a lot larger (c.f. Fig. in order to re- 
produce t 2 obs = 0.025 (i.e. (t 2 ) = 0.015). Although 
6 1 C Ori has shown evidence of variability in the stel- 
lar Hell absorption line (Conti 1972; Walborn 1981) 
and of a periodicity of 15 days in the Ha equivalent 
width (Stahl et al. 1993), no photometric variations 
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of the continuum on timescales from days to years 
have been reported (van Genderen et al. 1989). 

H II regions powered by Luminous Blue Variable 
(LBV) stars such as the Carina nebula are poten- 
tial candidates for finding a cause and effect rela- 
tion between t 2 and variability. One would expect 
i 2 to be significantly larger in nebulae excited by 
LBV stars. Unfortunately, we have no information 
at hand about t 2 in the Carina nebula for the pur- 
pose of comparison. 

A general and compelling argument against 
source variability as the main cause of the fluc- 
tuations is that for metal lines like Cm]AA1909, 
[Oiii]A5007, the bulk of their luminosity occurs at 
a nebular depth such that y < 10 -1 ' 7 , the fluctua- 
tions in temperature (due to variability) are negli- 
gible in that inner zone of the nebula despite sub- 
stantial changes in y, as shown in Fig. |2Jd. This 
means that (t 2 (0 +2 )) would be much smaller than 
(i 2 (H+)) calculated with model B" and A = 0.04. 
(Even larger values of A would lead to insignificant 
values of ^i 2 (0 +2 )).) In this particular model, it 
can be shown using mappings ic that 90% of [O in] 
flux is emitted in the inner regions where no fluctu- 
ations occur due to variability, as illustrated in Fig. 
03. Therefore our models, based on the variabil- 
ity of Qh, cannot realistically reproduce the large 
observed t 2 bs deduced from intermediate excitation 
emission lines. 

3.2. Variable shadows from moving opaque 
condensations 

There are other mechanisms than variable stars 
that can produce a transient variability in the prop- 
agating ionizing flux. For instance, the displacement 
of dense neutral condensations in planetary nebu- 
lae or H II regions (proplyds) can introduce tempo- 
ral variations in the ionization structure. Within the 
shadow, behind these opaque condensations, the gas 
is either neutral or partly ionized by the diffused field 
from the surrounding nebula. If these condensations 
of lateral size D had a tangential velocity component 
V (distinct from the nebular gas), the tangential dis- 
placement of the shadow inside the nebula would 
produce effects similar to that of an ionizing source, 
which is switched off, then on again, generating in 
turns an RF and IF, respectively. The important 
timescales in this problem are T rec and T cross , the 
crossover timescale (r c ,. oss = D/V) of the shadow 
over a distance given by the diameter of the opaque 
eclipsing condensation. We expect that the effect on 
the shadowed region would be strongest when T cross 
becomes comparable to r rec , since it would first in- 



TABLE 1 

FLUCTUATIONS FROM INTERMITENTLY 
SHADOWED IONIZED REGIONS 



n 


A a 


n x a 


(t 2 ) 


10 


A AO CT 

0.025 


0.25 


a A A o 

0.042 


5 


0.05 


0.25 


0.081 


2.5 


0.1 


0.25 


0.15 


1.25 


0.2 


0.25 


0.23 


10 


0.1 


1.0 


0.063 


10 


0.05 


0.5 


0.045 


10 


0.025 


0.25 


0.042 


10 


0.0125 


0.125 


0.021 


10 


0.00625 


0.0625 


0.013 



a Fraction of cycle during which the ionizing 
flux is turned off. 



duce a strong temperature drop and substantial re- 
combination followed by the propagation of hot IF 
after the "eclipse" has ended. On the global scale of 
the nebula, these shadowed regions would cause tem- 
perature fluctuations similar to those calculated for 
a star that is briefly turned off. Whether the tem- 
perature fluctuations would be felt over the whole 
nebula would depend on how many such shadowed 
regions fill the nebular volume. 

As an estimate of the t 2 expected in that case, we 
ran models in which the source was turned off dur- 
ing a time (II x A) T rec where II is the period of the 
off-on cycle in units of the recombination timescale 
and A is the fraction of the cycle during which the 
source is off. The idea behind using periodicity to 
approach this problem is that in this way we can 
crudely estimate the nebular fluctuations due to the 
shadows by associating A to the shadowed fraction, 
that is, to the effective covering factor of the source 
due to all the neutral condensations present. To a 
first order, this approach is validated by models for 
which we kept II x A = 0.25 constant, since we find 
that (i 2 ) cx A, at least when A is small (see first 4 
models in Table QJ. A sequence of models in which 
we change II x A is also given in Tabled Note that 
since the amplitude of the variations are extreme in 
the case of an on-off variability type, we expect that 
(i 2 (0 +2 )) would be substantial as well 11 . For the 

11 This is indirectly confirmed by the LINER model of Era- 
cleous, Livio & Binette (1995) who found a significantly higher 
[O III] A4363/[0 in]A5007 temperature sensitive ratio in their 
time-dependent model as a result of the very large variability 
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optimal regime where II x A ~ 1, a covering factor 
of 4% is required for the proplyds if we aim to reach 
t 2 obs = 0.025. Since these models were computed for 
a low value of U, we can estimate that (i 2 ) would 

be higher by factor about three, since (t 2 ) (x VU 
as discussed in Hence maybe a covering factor 
as low as 1% would suffice. However, the fraction 
of the Orion nebula volume shadowed by proplyds 
is estimated to be only 0.1% (William Henney: pri- 
vate communication), which therefore rules out this 
mechanism. 

4. CONCLUSIONS 

Our time-dependent models rule out the possibil- 
ity that the ionizing source variability be the cause 
of the variations because OB stars are not known 
to vary at the required 15% level (of Qh)- Further- 
more, even if they did, (t 2 (0+ 2 )) would still be much 
smaller than our computed (i 2 (H+)), making the ob- 
served values (c.f. Luridiana et al. 1999; Esteban et 
al. 1998) even more unattainable. As for the model 
of shadow crossings by opaque condensations, it re- 
quires a space density of proplyds a factor ten beyond 
that estimated by observations. 

An interesting and general result about cyclicly 
variable ionizing sources is that the nebula becomes 
"on average" less massive but somewhat hotter. De- 
spite the fluctuating size and temperature of the neb- 
ula, the time- averaged luminosity of all recombina- 
tion lines remains the same as for the steady-state 
case. This is also the case for the time-averaged en- 
ergy radiated through all the forbidden lines as com- 
pared to the steady-state case (individual forbidden 
line ratios must in general come out different since 
the nebula is hotter). 

The work of LB was supported by the CONA- 
CyT grant 32139-E. PF acknowledges support from 
the Region Rhone-Alpes. The comments received 
from William Henney about proplyds were particu- 
larly useful. We thank Jane Arthur for her contri- 
bution during the workshop on 3D-hydro in Mex- 
ico City (February 1999) during which the codes 
YGUANA and YGUAZU were written. We are also in- 
debted to AR who generously funded this event. 

A. yguana — THE EQUATIONS 

In this Appendix, we describe the set of equa- 
tions used to follow R-type ionization fronts within 
a spherical nebula photoionized by a variable source. 
The new code, hereafter called YGUANA, has been 

amplitudes of their source. 



written in Fortran 77 and will be made accessible to 
researchers upon request. 

We consider the problem of photoionization of a 
thick spherical gas shell by a central, time-varying 
source. We make the following simplifying assump- 
tions concerning the source and the gas distribution. 

• The ionizing radiation from the central source is 
taken to be monochromatic at a frequency z7 > 
vo (see value in Table 2) where hi'o= 13.6 eV is 
the ionization potential of hydrogen. The cen- 
tral source is considered point-like (i.e. its char- 
acteristic size is much smaller than the inner 
radius of the gas shell). 

• For the purpose of radiation transfer, the gas 
shell consists only of hydrogen. The gas density 
is static in time and its value is either a constant 
or a function of radius. 

• The atomic physics is simplified along the lines 
developped in IA.2l and IA.3I All essential phys- 
ical processes are considered (e.g. approximate 
cooling by metals), although the estimation of 
their rates is limited to first order approxima- 
tions. 

Given these assumptions, all variables of the 
problem depend only on the radius from the central 
source (spherical symmetry) and on the time. 

A.l. Time- dependent transfer equation 

The time-dependent equation of transfer for 
monochromatic ionizing radiation, and for spher- 
ically symmetric problem is (in spherical coordi- 
nates): 

1 dT 1 d(r 2 T) 

c at r z or 

where !F(r,t) is the ionizing photon flux (in photon 
s _1 cm -2 ) at a radius r from the central source (in 
cm) and at the time t (in s), n(r, t) is the local opacity 
at V (in photon s _1 cm -3 ), and c is the speed of 
light (in cm s _1 ). We have ignored the scattering of 
the ionizing radiation by dust and the generation of 
the diffuse field, which would have introduced 'local' 
source terms in equation ijlfl. 

For a pure hydrogen medium and monochromatic 
radiation, the opacity yields the relation: 

«(r,t) = a(V) en H (r) [1 - f(r,t)] F(r,t) (5) 

where a(u) is photoionization cross section of hydro- 
gen at 77 (in cm 2 ), e is the filling factor of the shell 
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(assumed uniform and static), nn(r) is the hydrogen 
number density (in cm -3 ), and f(r,t) is the ioniza- 
tion fraction of the hydrogen. The following relation 
has been used to compute a as a function of 77: 



77(77) = (,.:-! ■. 10 1N ( cm 2 



(6) 



where the value of v is chosen according to the type 
of ionizing source and is listed in Table 2. vq is the 
H ionizing threshold. A monochromatic treatment 
implies that the progressive hardening of the ioniz- 
ing radiation (and hence of the temperature) away 
from the central source cannot be taken into account. 
Furthermore, in the inner region of planetary nebu- 
lae and AGN, heating by photoionization of He + is 
substantial and lead to higher temperatures, this ef- 
fect cannot be taken into account in our simplified 
pure hydrogen nebula. 

A. 2. Time- dependent ionization balance of hydrogen 

To solve the transfer equation we need to 
know the ionization fraction f{r,t) of the hydrogen, 
and, therefore, to solve the ionization balance of hy- 
drogen. The ionization fraction yields the following 
equation: 



71(77) [1 - f(r,t)] T(r,t) 

V V ' 

photoionization 

+ 7 (T)n H (r) /(r,t) [1 - f(r,t)] 

collisional ionization 

- a B (T) n H (r) f 2 (r,t) 

S J ^ 

recombination 

(7) 

where T(r,t) is the temperature of the gas (in 
K), 7(T) is the collisional ionization coefficient (in 
s _1 cm 3 ), and «b(T) is the recombination rate to ex- 
cited states of hydrogen (in s _1 cm 3 ). We implicitly 
assume the on-the-spot approximation in our calcu- 
lations. The following analytical relations have been 
used to compute 7 and as a function of T: 



7(T) 



exp 



A + B x 



/ T 



V10 4 K 



— 1 3 

s cm 



(8) 



with A = — 19 and B = — 16, and 

/ rp \ -0.85 



«b(T) 



2.6 x 10 
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\10 4 K 



s^cm 3 (9) 



A. 3. Time- dependent energy balance 
The third major equation describing our system 
is the energy balance equation for the gas: 



^(r,*) = - 5 A(r,t) 



(10) 



TABLE 2 

PARAMETERS FOR THE RADIATIVE 
COOLING CURVE a AND THE VALUE OF V 



Environment 


A 






B 




hv (eV) 


Hll region 3.8 


x 10" 


24 


4.0 


x 10" 


24 


20 b 


Planetary nebula 4.0 


x 1Q- 


24 


3.5 


x 10" 


24 


40 c 


Active nucleus 4.2 


x 10" 


-24 


3.0 


x 10" 


24 


43 d 



a A and B are in units of erg s _1 cm 3 . 
b Stellar atmosphere of temperature 40 000 K. 
c Black body of temperature 150 000 K. 
d Power law of index —1.3 truncated at IkeV. 



where P is the gas pressure (in dyne cm 2 ) and 

A is the net cooling rate per unit of volume (in 

,-3 



erg s cm 3 ). A yields the relation: 



AM) 



£(T) 4(r) / 2 (r,t) 



radiative coolii 



+ tob 7 (T) ni(r)f(r,t) [1 - f(r,t)} 

S v ' 

cooling due to collisional ionization of H 

+ hv q(T) ni(r) f(r,t) [1 - f(r,t)] 

cooling due to collisional excitation of H 

- a(v) h(V-vo) n H (r) [1 - f(r,t)] T(r,t) 

heating due to photoionization of H 

(ii) 

where h is the Planck constant (6.63 x 10 27 erg s). 
The coefficients £(T), 7 (T), q(T) and 75(77) are given 
in Eq. ©, JEl and ©, respectively. 

In order to have a realistic energy balance for the 
gas, the radiative cooling term includes the losses 
due to metals. £(T) has been computed using the 
following relation, which is an analytical fit to the 
cooling function obtained for a photoionized gas us- 
ing MAPPINGS ic and assuming solar metallicities: 



£(T) = A 



B 



T 



10 4 K 



2.5 



erg s 



_1 cm 3 



(12) 



where the values of A, B depend on the type of envi- 
ronment. Three sets of values are listed in Table 2. 

This radiative cooling term applies only to a fully 
ionized medium. We do not take into account how 
this term varies with the ionization parameter al- 
though we consider separately the cooling due to ex- 
citation or ionization of H°, which can be significant 
across ionization fronts. To compute the net cooling 
rate we therefore added to Eq. 1|11|) the cooling terms 
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due to both collisional ionization and collisional ex- 
citation of hydrogen as well as the heating term due 
to photoionization of hydrogen. The collisional exci- 
tation rate q(T) yields the following relation [similar 
to Eq. ©]: 



3(T) 



exp 



C + D x 



s ^m 3 



10 4 K, 

(13) 

with C = — 18 and D = —11. The determination 
of the constants A, B, C and D are based on the 
coefficients found in Osterbrock (1989). 

A. 4. Equation of state 

In order to close our system of equations, we need 
to express the equation of state of the gas. We have 
used the perfect gas equation of state: 



P(r,t) = [l + /(r,t)] n H (r)fcT(r,t) 



(14) 



where k is the Boltzmann's constant (1.38 x 10~ 16 
erg K" 1 ). From Eq. 1)1(1 1) and H14fl . we can derive the 
equation for the temperature: 



1 



fcn H (r) [1 + /M)j 



-§ A aet (r,t) - n H (r) A; T(r,t) |f (r,t) 

'(15) 

Together, equations J7J and (|15|) form a 

closed system of equations. 

B. yguana — THE ALGORITHM 

In this Appendix, we describe the practical im- 
plementation of a code to solve the set of equations 
describe inlAl 

B.l. Boundary conditions 

To solve the problem, we need to define a set 
of boundary conditions. First, we have to assume 
initial (i.e. at t = s) temperature and ionization 
structures for the whole shell (i.e. from its inner ra- 
dius, Ri n , to its outer radius, R ou t), as well as the 
corresponding, initial distribution of ionizing radia- 
tion: 



Vr, R in < r < R 



/(r,0) - fair) 
T(r,0) = T (r) (16) 



Last, we need to know the ionizing photon flux 
of the source Source at R in as a function of time: 



Vt, jF(Ri n ,i) — !F BOUICe (t) 



(17) 



Assuming that the cavity inside the shell (r < Ri n ) 
is empty, ^source (t) can be expressed as a function of 
the central source ionizing photons production rate, 
Qsourcc (in photon s _1 ): 



Qs 



4tt R 2 



(18) 



B.2. Algorithm 

In the following text, the indices i and j refer to 
the spatial and time axes, respectively. 

B.2.1. Transfer equation 

Once linearized, equation (J1J allows us to com- 
pute the ionizing photon flux Ti j at a radius r*j and 
time tj, as a function of the ionization and temper- 
ature structures of the shell and of the ionizing pho- 
ton flux at radii r^_i and r^, and time tj—\- A sketch 
of the progression of the computation in the spatial 
and temporal grids is shown in Fig. We have the 
following recursive equation for J-'ij: 



CSt -77 



x [l-ti-ij-iSr-^ 



incoming photon flux (diluted and absorbed) 



i-1,3-1 



out-coming photon flux 

a(F)en H (r l _i) (1 - /.,-ij-i) 



where, 5r = fj — rj— i, 5t — tj — tj—i, and c is the 
speed of light. For this equation to be valid, we must 
always have: 



Sr 



< 1 



*■ « (— ) cVi 



(20) 



If 5t is larger than the typical atomic physic time 
scales (e.g. the recombination time), the timestep 
is divided in smaller time intervals. This prevents 
(i.e. the opacity) from changing significantly 
during the timestep. 

B.2.2. Ionization balance 

Once J-'ij is computed, the time-dependent ion- 
ization balance [equation Q] is then solved. The 
new ionization fraction at tj, fij, is computed as a 
function of fij-i, Tjj_i and (JFij + J 7 i.j-i)/2. The 
algorithm used to solve equation (JJJ is the same as 
in MAPPINGS ic (see Appendix in Binette, Dopita & 
Tuohy 1985). 
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0-1 ,7-1) 



8r 



r - r 



O'-i ,j) 



0,7-1) 




CO 



,7) 



Fig. 6. Sketch of the progression of the computation 
along the spatial (i index) and temporal (j index) grids. 
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B.2.3. Temperature equation 

At last, the temperature Tij is computed us- 
ing the following recursive equation [as derived from 
Eq. JEJ]: 

Tjj-i — StQi t j—i — 

3kn H (ri) (l+fi.j) ( 21 ) 

0jj and represent the change in temperature 
due to the net cooling and the change in 'molecu- 
lar' weight, respectively. If necessary, the ionization 
balance and the temperature equation were solved 
iteratively. 
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